# Relative error estimates for the initial state in %
error_stats = tibble::tibble(
sV = 0.26, # Glacier volume
sA = 0.13, # Glacier area
sTh = 0.13, # Glacier thickness
sMelt = 2) # Glacier melt 11 Modeling of Discharge from Glacier Melt
The glacier melt modeling in RSMinerve (Garcia Hernandez et al. 2020) is done (at the time of writing, March 2022) with the GSM model which features a constant area and an unlimited glacier reservoir. It is suitable for short-term simulations of glacier melt where effects of glacier volume change can be neglected but it is not well suited for climate change impact studies where substantial changes in glacier volume are to be expected. However, by assuming that discharge from snow melt can be treated independently from discharge from glacier melt, the discharge contribution from glacier melt can be simulated in R and included as source into RSMinerve models. The following chapter relies on the basic understanding of the glacier mass balance as presented and discussed in Chapter Snow and Glacier Data, describes the glacier modelling tools in the package riversCentralAsia and demonstrates how to use them for joint applied hydrological modelling with RSMinerve.
11.1 Temperature Index Model
One of the arguably simplest models to describe glacier melt (or snow melt) is the temperature index model. Hock (2003) describes several variants of the temperature index model for simulating glacier melt. The riversCentralAsia package implements the simplest temperature index melt model described in (Hock 2003) in the function glacierMelt_TI (Equation 11.1), requiring only temperature time series and 2 parameters as input.
\[ M = \biggl\{ \begin{array}{l, l} 0, & T < T_{threshold} \\ f_{M} \cdot \left( T - T_{threshold} \right), & T >= T_{threshold} \end{array} \qquad(11.1)\]
where \(M\) is the glacier melt in \(mm/d\), \(T\) is the daily average temperature in \(^{\circ} C\). The two parameters \(f_{M}\) and \(T_{\text{threshold}}\) refer to the melt factor and the threshold temperature above which glacier melt occurs and need to be calibrated. They have the units \(\frac{mm}{^{\circ} C \cdot d}\) and \(^{\circ} C\) respectively. Glacier melt is calculated in daily time steps.
11.2 Glacier Mass Balance
The glacier mass balance is simplified to Equation 11.2:
\[ \Delta S = P-M = A_{\text{imbal}} \qquad(11.2)\]
where the change of water storage (\(\Delta S\)) is equal to the precipitation (\(P\)) minus the glacier melt (\(M\)). Typically melt exceeds precipitation and we have negative \(\Delta S\), that is imbalance ablation, indicating glacier storage loss. The glacier mass balance is calculated in annual time steps. We thereby refer to the hydrological year starting on October 1st of the previous year to take a full accumulation and ablation season into account.
11.3 Glacier Volume Development
As glaciers melt, their volume changes. This has to be taken into account for the long-term simulation of glacier discharge. To determine the initial glacier volume, the area of the geometry of the Randolph Glacier Inventory (RGI) v6.0 data set is multiplied with the average thickness of the glacier (the Farinotti data set). From the combination of these two data sets, the well established area-volume relationship by Erasov (1968) can be verified (see section on Area-Volume scaling).
Please note that the data and the retrieval of the data are described in Part II of this book. Glaciers larger than 1km2 are sub-divided into elevation bands of 100 m altitude to account for elevation dependent temperature forcing. The large glaciers melt from the lowest elevation band to the highest elevation band whereby the glacier melt is subtracted from the glacier volume of lowest elevation band that is still glacierized. The small glaciers are not spatially discretized and thus they are melted homogeneously.
For annual time step \(t\), the evolution of the glacier volume is calculated as follows:
\[ A(t) = \text{glacierArea RGIF}\bigl(V(t)\bigr) \qquad(11.3)\]
\[ Q_{glacier}(t) = M(t) \cdot A(t) \qquad(11.4)\]
\[ V(t+1) = V(t) + \Delta S = V(t) + \text{glacierImbalAbl}\bigl(M(t)\bigr) \qquad(11.5)\]
\(Q_{\text{glacier}}\) can be calibrated against glacier discharge derived from the Miles & Hugonnet data sets. The automated calibration is currently not included in the riversCentralAsia package.
The function glacierArea_RGIF() is an empirical scaling function analogue the inverse of the scaling function derived by Erasov (1968) but based on the modern RGI v6.0 glacier geometries and the glacier thickness data set by Farinotti et al. (2019). Section Area-volume scaling shows the derivation of the scaling relationship and its comparison to the relationship proposed by Erasov (1968). The package riversCentralAsia implements volume-area and area-volume scaling functions based on both, Erashov and RGI-Farinotti data, allowing the estimation of glacier areas based on glacier volumes (glacierArea_Erasov and glacierArea_RGIF) and estimations of glacier volumes based on glacier areas (glacierVolume_Erasov and glacierVolume_RGIF).
The function glacierImbalAbl is an empirical scaling function relating glacier imbalance ablation, i.e. the glacier melt from storage loss. It is derived from the glacier discharge data set by Miles et al. (2021) and the glacier thinning data set by Hugonnet et al. (2021) (see section on glacier discharge-thinning scaling for the derivation of the relationship).
11.4 Summary of Workflow
To model the contribution to discharge from glacier melt in a basin, the following steps are required:
- Pre-processing of GIS layers
- Pre-processing of climate forcing data
- Calculation of daily glacier melt
- Calculation of annual glacier masss balance
- Scaling of annual glacier contribution to daily values
- Aggregation of per glacier contribution to sub-basins (optional)
- Writing of RSMinerve source intput files
- Integration of glacier discharge sources in RSMinerve
11.5 Propagation of Uncertainty
The initial glacier volume of each glacier is attributed an uncertainty of p/m 26%. This number is based on the average uncertainty of the glacier volume per RGI region reported in Farinotti et al. (2019). The uncertainty of the area of the RGI v6.0 glacier outlines is not known. It is therefore assumed that the error of the glacier volume stems to 50% from the estimation of the glacier thickness and to 50% from the glacier area. We further assume that the errors of the glacier area and glacier thickness are un-correlated and can thus estimate the uncertainties of the glacier area data and the glacier thickness data to be p/m 13% each.
For the non-linear relationships in glacierVolume_RGIF and glacierArea_RGIF, the standard deviation of the residuals of the fit was computed. The estimated error of the fit is assumed to be equal plus/minus twice the standard deviation of the residuals and yields 31% and 53% respectively. The residuals are not normally distributed and their actual distribution is unknown. Further, uncorrelated errors and the applicability of linear error propagation are assumed. Therefore, the error of the function outputs is simply computed by adding the error of the function input to the error of the fit.
\[ \varepsilon_{V} = \varepsilon_{A} + \varepsilon_{\text{glacierVolume RGIF}} = 0.26 + 0.31 = 0.57 \qquad(11.6)\]
\[ \varepsilon_{A} = \varepsilon_{V} + \varepsilon_{\text{glacierArea RGIF}} = 0.26 + 0.53 = 0.79 \qquad(11.7)\]
Error estimates for the temperature index model are not available. A conservative relative error of 2 is therefore assumed, indicating that the estimated glacier melt is within a range of plus/minus 2 times it’s value.
The error of the imbalance ablation amounts to 73%. Assuming independent errors from the temperature index model for glacier melt and the scaling function between imbalance ablation and glacier melt, the error of the imbalance ablation amounts to approximately:
\[ \varepsilon_{Q_{\text{imb,melt}}} = \varepsilon_{Q_{\text{M}}} + \varepsilon_{\text{glacierImbalAbl}} = 2 + 0.73 = 2.73 \qquad(11.8)\]
The errors stated above relate to initial estimates of the glacier areas, volumes, total ablation and imbalance ablation. The errors of the glacier model variables of each subsequent time step depend on the errors of the previous time steps, i.e. they are not uncorrelated over time. For simplicity reason, this non-linear propagation of errors is neglected. In any case, the uncertainty margins for any glacier melt modelling done with the presented method are large. The following table summarises the estimated errors for each variable
11.6 Area-volume scaling
library(sf, quietly = TRUE)
library(tidyverse, quietly = TRUE)
library(gridExtra, quietly = TRUE)
library(broom, quietly = TRUE)
devtools::install_github("hydrosolutions/riversCentralAsia", quiet = TRUE)
library(riversCentralAsia, quietly = TRUE)11.6.1 Load the data
# Load the RGIv6.0 data set of the RGI region of Central Asia
rgi <- st_read("../caham_data/central_asia_domain/glaciers/RGIv60/13_rgi60_CentralAsia.shp",
quiet = TRUE) |>
sf::st_make_valid()
rgi <- rgi |>
dplyr::select(RGIId) |>
mutate(Area_m2 = as.numeric(st_area(rgi)))
# Load the glacier thickness and extract avearage glacier thickness per glacier
load("../caham_data/central_asia_domain/glaciers/Farinotti/glacier_thickness_extractedRGIreg13.RData")
temp <- as.numeric(st_area(rgi))
rgi <- rgi |>
left_join(glacier_thickness |>
dplyr::select(RGIId, mean), by = "RGIId") |>
rename(thickness_m = mean) |>
mutate(Volume_m3 = Area_m2 * thickness_m)
# Prepare the data for curve fitting. Exclude Fedchenko glacier as it is
# so much larger than all other glaciers in RGI region 13.
data <- rgi |>
st_drop_geometry() |>
dplyr::select(RGIId, Area_m2, thickness_m, Volume_m3) |>
mutate(Area_km2 = Area_m2 * 10^(-6),
Volume_km3 = Volume_m3 * 10^(-9)) |>
dplyr::filter(Area_km2 < 300) |>
drop_na()
# Fit a model of similar shape as Erasovs area-volume relationship.
# V = 0.027 * A ^1.5
# Note: We use the expected error of the glacier thickness estimate of 0.26
# (as reported by Farinotti et al., 2019) to weight the least squares.
err <- 0.26
nlVAw <- nls(Volume_km3 ~ a * Area_km2^b,
data = data,
weights = data$Area_km2*0+1/err^2,
start = list(a = 0.05, b = 1.2))
# Add the modeled volume to the "observed" volume.
data_nlVAw <- data |>
mutate(.fitted = augment(nlVAw)$.fitted,
.fitted_km3 = .fitted,
.resid = augment(nlVAw)$.resid,
.relresid = .resid/.fitted_km3,
.se = 2*sd(.relresid),
.lb = ifelse(.fitted_km3*(1- .se) < 0, 0, .fitted_km3*(1- .se)),
.ub = .fitted_km3*(1+ .se),
Volume_Erasov_km3 = glacierVolume_Erasov(Area_km2))
# Plotting the results
a <- ggplot(data_nlVAw) +
geom_ribbon(aes(Area_km2, .fitted_km3, ymin = .lb, ymax = .ub),
colour = NA, alpha = 0.2, fill = "forestgreen") +
geom_point(aes(Area_km2, Volume_km3, colour = "Obs"), size = 0.4,
alpha = 0.4) +
geom_point(aes(Area_km2, Volume_Erasov_km3, colour = "Erasov"), size = 0.4,
alpha = 0.4) +
geom_point(aes(Area_km2, .fitted_km3, colour = "Sim"),
alpha = 0.4, size = 0.4) +
scale_colour_manual(name = "Legend",
values = c("Obs" = "black", "Erasov" = "green",
"Sim" = "forestgreen")) +
ylab(expression("Volume [" * km^3 * "]")) +
xlab(expression("Area [" * km^2 * "]")) +
theme_bw()
b <- ggplot(data_nlVAw) +
geom_ribbon(aes(Area_km2, .fitted_km3, ymin = .lb, ymax = .ub),
colour = NA, alpha = 0.2, fill = "forestgreen") +
geom_point(aes(Area_km2, Volume_km3, colour = "Obs"), size = 0.4,
alpha = 0.4) +
geom_point(aes(Area_km2, Volume_Erasov_km3, colour = "Erasov"), size = 0.4,
alpha = 0.4) +
geom_point(aes(Area_km2, .fitted_km3, colour = "Sim"),
alpha = 0.4, size = 0.4) +
scale_colour_manual(name = "Legend",
values = c("Obs" = "black", "Erasov" = "green",
"Sim" = "forestgreen")) +
ylab(expression("Volume [" * km^3 * "]")) +
xlab(expression("Area [" * km^2 * "]")) +
ylim(c(0, 4)) +
xlim(c(0, 20)) +
theme_bw()
grid.arrange(a, b, nrow = 1)Warning: Removed 236 rows containing missing values (geom_point).
Removed 236 rows containing missing values (geom_point).
Removed 236 rows containing missing values (geom_point).
For small glaciers with an area \(< 10 km^2\) Erasov (\(V_{\text{Erasov}} = 0.027 \cdot A^{1.5}\)) performs better than the newly fitted data from RGI and Farinotti (\(V_{\text{RGIF}}=\) 0.061 \(\cdot A^{b}\) with \(b=\) 1.142). For glaciers larger than \(10 km^2\) area, the newly fitted function glacierVolume_RGIF should be used.
11.7 Glacier discharge-thinning scaling
11.8 Demonstration
We demonstrate the above described method with the data from the Atbashy basin.
Warning
All required data is available from the downloadable data package. Note however that this folder required 12 GB of local storage space!
library(tmap, quietly = TRUE)
library(sf, quietly = TRUE)
library(raster, quietly = TRUE)
library(exactextractr, quietly = TRUE)
library(tidyverse, quietly = TRUE)
library(lubridate, quietly = TRUE)
devtools::install_github("hydrosolutions/riversCentralAsia", quiet = TRUE)
library(riversCentralAsia, quietly = TRUE)
# Path to the data directory downloaded from the download link provided above.
# Here the data is extracted to a folder called atbashy_glacier_demo_data
data_path <- "../caham_data/SyrDarya/Atbashy/"11.8.1 Forcing
There will be a separate Section to demonstrate how to prepare the forcing. For now, we provide you with the pre-processed forcing data and give you just a brief overview over the data source.
Historical Temperatures
As meteorological data for high elevations in Central Asia is very scarce we use the CHELSA v2.1 data set (Karger et al. 2017). This is a global data set of forcing data for hydrolgical models based on ERA5 but corrected for biases and especially suited for high elevations. The daily CHELSA forcing has been cut to the Central Region by the originator of the data set, D. Karger, WSL and extracted to the hydrological response units of the glaciers in the Atbashy basin by T. Siegfried.
Future Temperatures
Future temperature development per glacier or elevation band is extracted from the 3 CMIP6 GCM models with highest priorities for the region downloaded from COPERNICUS. We take 4 socioeconomic scenarios into account, covering 4 different emission scenarios. The temperatures of the climate models are bias corrected using the CHELSA data and a quantile mapping method. More details in the climate data preparation section.
fut_sim <- readRDS(file = paste0(data_path, "CLIMATE/fut_sim_glacier_tas_qmapped.RDS"))
# Plot the temperature time series for a given glacier/elevation band
glacier <- "RGI60-13.08930_1"
# Extract the temperature for the selected glaciers for all GCMs and SSPs
scenarios <- names(fut_sim)
fut_temp <- NULL
for (scenario in scenarios) {
fut_temp <- rbind(fut_temp,
fut_sim[[scenario]] |>
dplyr::select(Date, all_of(glacier)) |>
mutate(Scenario = scenario))
}
fut_temp <- fut_temp |>
mutate(Hyear = hyear(Date)) |>
group_by(Hyear, Scenario) |>
summarise(Date = first(Date),
Temp = mean(get(glacier))) |>
separate(Scenario, into = c("GCM", "SSP"), sep = "_") |>
ungroup() |>
dplyr::filter(Hyear > min(Hyear) & Hyear < max(Hyear),
GCM != "IPSL-CM6A-LR")
# Plot annual data for better readability
ggplot() +
geom_line(data = hist_obs |>
mutate(Hyear = hyear(date)) |>
group_by(Hyear) |>
summarise(date = first(date),
Temp = mean(get(glacier))) |>
ungroup() |>
dplyr::filter(Hyear > min(Hyear) & Hyear < max(Hyear)),
aes(date, Temp)) +
geom_line(data = fut_temp, aes(Date, Temp, colour = GCM,
linetype = SSP)) +
scale_colour_viridis_d() +
labs(x = "Date", y = "T [deg C]") +
theme_bw()
11.8.2 Glacier Geometry
# Load data
## Digital elevation model (DEM)
dem <- raster(paste0(data_path, "GIS/16076_DEM.tif"))
## Glacier outlines from the Randolph Glacier Inventory (RGI) v6.0
rgi <- st_read(paste0(data_path, "GIS/16076_Glaciers_per_subbasin.shp"),
quiet = TRUE) |>
st_transform(crs = crs(dem))
## Outlines of hydrological response units for the modelling of glacier discharge.
rgi_elbands <- st_read(paste0(data_path,
"GIS/rgi_glaciers_atbaschy_el_bands.shp"),
quiet = TRUE) |>
st_transform(crs = crs(dem))
# Load a pre-processed raster file with glacier thickness in the Atbashy basin
# (see vignette [glaciers 01]{glaciers-01-intro.html} for details on the pre
# -processing)
glacier_thickness <- raster(
paste0(data_path,
"GLACIERS/Farinotti/pre-processed_glacier_thickness.tif"))
# Load the glacier thickness data set and filter it to the glaciers in the
# catchment of interest.
hugonnet <- read_csv(paste0(
data_path, "/GLACIERS/Hugonnet/dh_13_rgi60_pergla_rates.csv"))
hugonnet <- hugonnet |>
dplyr::filter(rgiid %in% rgi$RGIId) |>
tidyr::separate(period, c("start", "end"), sep = "_") |>
mutate(start = as_date(start, format = "%Y-%m-%d"),
end = as_date(end, format = "%Y-%m-%d"),
period = round(as.numeric(end - start, units = "days")/366))
glaciers_hugonnet <- rgi |>
left_join(hugonnet |> dplyr::select(rgiid, area, start, end, dhdt, err_dhdt,
dvoldt, err_dvoldt, dmdt, err_dmdt,
dmdtda, err_dmdtda, period),
by = c("RGIId" = "rgiid"))
# Only keep the variables we need for this analysis
rgi_elbands <- rgi_elbands |>
dplyr::select(RGIId, Area, elvtn_b) |>
rename(Area_tot_glacier_km2 = Area) |>
mutate(ID = paste0(RGIId, "_", elvtn_b))
# Get mean elevation of each glacier/elevation band from DEM
rgi_elbands$z_masl <- exact_extract(dem, rgi_elbands, "mean", progress = FALSE)
# Update the glacier area within the basin boundaries
rgi_elbands$A_km2 <- as.numeric(st_area(rgi_elbands))*10^(-6)
glaciers <- unique(rgi_elbands$RGIId)
# Calculate the average glacier thickness for each elevation band.
rgi_elbands$thickness_m = exact_extract(glacier_thickness,
rgi_elbands, "mean", progress = FALSE)11.8.3 Processing of Forcing
Temperature time series for each glacier/elevation band are available. Figure 11.3} shows temperature time series extracted from the CHELSA data set on the example of one of the larger glaciers in the Atbashy basin. The data is aggregated to the hydrological year to better illustrate the temperature gradient over the elevation. The highest temperatures are measured at the lowest elevations and vice versa.
ggplot(hist_obs[, 1:9] |>
pivot_longer(-date, names_to = "ID", values_to = "Temp") |>
mutate(Hyear = hyear(date)) |>
group_by(Hyear, ID) |>
summarise(date = first(date),
Temp = mean(Temp))|>
ungroup() |>
dplyr::filter(Hyear>min(Hyear) & Hyear<max(Hyear)) |>
left_join(rgi_elbands |>
st_drop_geometry() |>
dplyr::select(ID, z_masl), by = "ID") |>
mutate("Elevation [masl]" = as.factor(round(z_masl)))) +
geom_line(aes(date, Temp, colour = `Elevation [masl]`)) +
scale_colour_viridis_d() +
labs(x = "Date", y = "T [deg C]") +
theme_bw()
11.8.4 Calculating Glacier Melt
The code snippet below illustrates how to calculate the annual glacier melt of the catchment. Figure 11.4 shows daily melt rates for the different elevation bands of one of the larger glaciers (same as in the figure above). The highest melt occurs at low elevation bands.
MF_small = 1
MF_large = 0.5
threshold_temperature = 0
Area <- rgi_elbands |>
st_drop_geometry() |>
dplyr::select(ID, A_km2) |>
pivot_wider(names_from = ID, values_from = A_km2)
# Assign different melt factors to large and small glaciers.
MF <- Area |>
mutate(across(everything(), ~MF_large),
across(ends_with("_1"), ~MF_small))
melt <- glacierMelt_TI(temperature = hist_obs |> dplyr::select(-date),
MF = MF,
threshold_temperature = threshold_temperature)
melt <- as_tibble(melt) |>
mutate(date = hist_obs$date) |>
relocate(date, .before = where(is.numeric))
ggplot(melt[, 1:9] |>
pivot_longer(-date, names_to = "ID", values_to = "Melt") |>
left_join(rgi_elbands |>
st_drop_geometry() |>
dplyr::select(ID, z_masl), by = "ID") |>
mutate("Elevation [masl]" = as.factor(round(z_masl)))) +
geom_line(aes(date, Melt, colour = `Elevation [masl]`)) +
scale_colour_viridis_d() +
labs(x = "Date", y = "Melt [mm/d]") +
theme_bw()
11.8.5 Compare to Measured Glacier Melt
Aggregate the daily glacier melt to annual data and compare it to the observed glacier melt. Note that the glacier melt simulated with the temperature index model (sim) is never negative whereas the glacier mass change derived from Hugonnet et al. (2021) and Miles et al. (2021) can be negative, indicating glacier growth.
melt_a_eb <- melt |>
pivot_longer(-date, names_to = "ID", values_to = "Melt") |>
mutate(Hyear = hyear(date)) |>
group_by(Hyear, ID) |>
summarise(date = first(date),
Melt = sum(Melt),
.lb_Melt = ifelse(Melt*(1-error_stats$sMelt)<0, 0,
Melt*(1-error_stats$sMelt)),
.ub_Melt = Melt*(1+error_stats$sMelt))|>
ungroup()
melt_a <- melt_a_eb |>
separate(ID, into = c("RGIId", "elB"), sep = "_") |>
group_by(Hyear, RGIId) |>
summarise(date = as_date(first(date)),
Melt = sum(Melt),
.lb_Melt = ifelse(Melt*(1-error_stats$sMelt)<0, 0,
Melt*(1-error_stats$sMelt)),
.ub_Melt = Melt*(1+error_stats$sMelt)) |>
ungroup() |>
dplyr::filter(Hyear > min(Hyear) & Hyear < max(Hyear))
glaciers_hugonnet <- glaciers_hugonnet |>
mutate(Qgl_m3a = glacierDischarge_HM(dhdt),
.lb_Qgl_m3a = ifelse(
dhdt > 0,
ifelse(Qgl_m3a*(1-error_stats$sQglgrowth)<0, 0,
Qgl_m3a*(1-error_stats$sQglgrowth)),
ifelse(Qgl_m3a*(1-error_stats$sQglmelt)<0, 0,
Qgl_m3a*(1-error_stats$sQglmelt))),
.ub_Qgl_m3a = ifelse(dhdt > 0,
Qgl_m3a*(1+error_stats$sQglgrowth),
Qgl_m3a*(1+error_stats$sQglmelt)))
melt_obs_a <- glaciers_hugonnet |>
dplyr::filter(RGIId %in% glaciers[6:9],
period == 1)
ggplot() +
geom_ribbon(data = melt_a |>
dplyr::filter(RGIId %in% glaciers[6:9]),
aes(date, Melt/1000, ymin = .lb_Melt/1000, ymax = .ub_Melt/1000,
fill = RGIId), colour = NA, alpha = 0.2) +
geom_ribbon(data = melt_obs_a |>
dplyr::filter(RGIId %in% glaciers[6:9]),
aes(start, -dmdtda,
ymin = -dmdtda-err_dmdtda,
ymax = -dmdtda+err_dmdtda,
colour = RGIId, linetype = "obs", fill = RGIId),
size = 0.2, alpha = 0.2) +
geom_line(data = melt_a |>
dplyr::filter(RGIId %in% glaciers[6:9]),
aes(date, Melt/1000, colour = RGIId, linetype = "sim")) +
geom_line(data = melt_obs_a, aes(start, -dmdtda, colour = RGIId,
linetype = "obs")) +
labs(x = "Date", y = "Melt [m weq/a]") +
scale_linetype_manual(name = "Source",
values = c("sim" = 1, "obs" = 2)) +
scale_colour_viridis_d() +
scale_fill_viridis_d() +
ggtitle(paste0("MF: ", MF, "Tth: ", threshold_temperature)) +
theme_bw()
The temperature index model can be calibrated with the specific glacier volume change provided by Hugonnet et al. (2021) (see Figure 11.5}). For the moment, manual calibration of the parameters is required.
11.8.6 Glacier Mass Balance
The following figure shows the components of the glacier mass balance for a few glaciers in the Atbashy basin.
glacier_balance <- glacierBalance(melt_a_eb = melt_a_eb,
rgi_elbands = rgi_elbands)
glacier_balance <- glacier_balance |>
mutate(.lb = ifelse(Variable == "A_km2",
ifelse(Value*(1-error_stats$sA)>0,
Value*(1-error_stats$sA), 0),
ifelse(Variable == "V_km3",
ifelse(Value*(1-error_stats$sV)>0,
Value*(1-error_stats$sV), 0),
ifelse(Variable == "Q_m3a",
ifelse(Value*(1-error_stats$sQglmelt)>0,
Value*(1-error_stats$sQglmelt), 0),
ifelse(Variable == "Qimb_m3a",
Value*(1-error_stats$sImbal), NA)))),
.ub = ifelse(Variable == "A_km2",
Value*(1+error_stats$sA),
ifelse(Variable == "V_km3",
Value*(1+error_stats$sV),
ifelse(Variable == "Q_m3a",
Value*(1+error_stats$sQglmelt),
ifelse(Variable == "Qimb_m3a",
Value*(1+error_stats$sImbal), NA)))))
ggplot(glacier_balance |>
dplyr::filter(RGIId %in% glaciers[6:9],
Hyear > min(Hyear) & Hyear < max(Hyear),
Variable %in% c("A_km2", "V_km3", "Q_m3a", "Qimb_m3a"))) +
geom_ribbon(aes(Hyear, ymin = .lb, ymax = .ub, fill = RGIId),
alpha = 0.2, colour = NA) +
geom_line(aes(Hyear, Value, colour = RGIId)) +
facet_wrap("Variable", scales = "free_y") +
scale_colour_viridis_d() +
scale_fill_viridis_d() +
theme_bw()
We now have glacier discharge (Q_m3s) and the unsustainable contribution to glacier discharge, the imbalance ablation (Qimb_m3a) which is negative for glacier loss and positive for growing glaciers. We are only interested in the contribution of imbalance ablation to river discharge, that is, only the negative part of Qimb_m3a is relevant to us.
11.9 From Annual to Daily Melt
The glacier mass balance is done on a yearly basis (if not at lower frequency). Hydrological models, however, typically run at higher frequency, for example monthly or daily time steps. One simple method to distribute glacier discharge on a year is to scale it according to the daily melt computed.
# Aggregate the daily melt per glacier
melt_mmd <- melt |>
pivot_longer(-date, names_to = "ID", values_to = "melt") |>
separate(ID, into = c("RGIId", "elB"), sep = "_") |>
group_by(date, RGIId) |>
summarize(melt = sum(melt)) |>
ungroup() |>
rename(M_mmd = melt)
# Compute the annual melt per glacier
melt_mma <- melt_mmd |>
mutate(Hyear = hyear(date)) |>
group_by(Hyear, RGIId) |>
summarise(M_mma = sum(M_mmd)) |>
ungroup()
# Rescale the annual imbalance glacier ablation with the daily melt rates
imbalAbl_m3s <- melt_mmd |>
mutate(Hyear = hyear(date)) |>
left_join(melt_mma, by = c("RGIId", "Hyear")) |>
left_join(glacier_balance |>
dplyr::filter(Variable == "Qimb_m3a") |>
transmute(Hyear = Hyear,
RGIId = RGIId,
Qimba_m3s = -1* Value/(60*60*24*365)) |>
mutate(Qimba_m3s = ifelse(Qimba_m3s < 0, 0, Qimba_m3s)),
by = c("RGIId", "Hyear")) |>
mutate(Qimb_m3s = Qimba_m3s * (M_mmd/M_mma),
Qimb_m3s = ifelse(is.na(Qimb_m3s), 0, Qimb_m3s))
# Visualise daily imbalance ablation
ggplot(imbalAbl_m3s |>
dplyr::filter(RGIId %in% unique(glacier_balance$RGIId)[1:7])) +
geom_line(aes(date, Qimb_m3s, colour = RGIId)) +
scale_colour_viridis_d() +
labs(x = "Date", y = "Glacier discharge from imbalance ablation [m3/s]") +
theme_bw()
12 Aggregation to Sub-Basins
Glacier discharge can be aggregated to sub-basin in RSMinerve. For this, a GIS layer needs to be prepared manually in QGIS because typically, the RGI v6.0 glacier outlines are not consistent with the river basin boundaries derived from SRTM. In QGIS the sub-basin layer is intersected with the RGI layer and the boundaries between the sub-basins are manually adjusted to be consistent with the glacier outlines. This step is important as any glaciers which are cut by sub-basin boundaries and thus are present in two different sub-basins will be counted doubly.
dem <- raster(paste0(data_path, "GIS/16076_DEM.tif"))
basin <- st_read(paste0(data_path, "GIS/16076_Basin_outline.shp"), quiet = TRUE)
hru <- st_read(paste0(data_path, "GIS/16076_HRU.shp"), quiet = TRUE) |>
st_make_valid() |>
mutate("Sub-basin" = gsub("_\\d*", "", name),
"Sub-basin" = gsub("Subbasin", "", `Sub-basin`))
rgi <- st_read(paste0(data_path, "GIS/16076_Glaciers_per_subbasin.shp"),
quiet = TRUE) |>
st_transform(crs = crs(dem))
# Visualise glaciers coloured by sub-basin
tmap_mode("view")
tmap_options(check.and.fix = TRUE)
#tm_shape(dem) +
# tm_raster(n = 6,
# palette = terrain.colors(6),
# alpha = 0.8,
# legend.show = TRUE,
# title = "Elevation (masl)") +
tm_shape(hru |> st_make_valid()) +
tm_polygons(col = "Sub-basin", lwd = 0.6) +
tm_shape(basin) +
tm_borders(col = "black", lwd = 0.6) +
tm_shape(rgi |> mutate("Glaciers" = gsub("_Subbasin", "", name_2))) +
tm_polygons(col = "Glaciers", lwd = 1, border.col = "black")The code junk below shows how to aggregate the daily per-glacier imbalance ablation rates to sub-basins.
Qimb_m3s_sub <- imbalAbl_m3s |>
dplyr::select(RGIId, date, Qimb_m3s) |>
left_join(rgi |>
st_drop_geometry() |>
dplyr::select(RGIId, name_2)) |>
group_by(date, name_2) |>
summarise(Qimb_m3s = sum(Qimb_m3s, na.rm = TRUE)) |>
ungroup()
ggplot(Qimb_m3s_sub |>
mutate("Sub-basin" = gsub("_Subbasin", "", name_2))) +
geom_line(aes(date, Qimb_m3s, colour = `Sub-basin`), alpha = 0.8) +
scale_colour_viridis_d() +
labs(x = "Date", y = "Glacier discharge from imbalance ablation [m3/s]") +
theme_bw()
12.1 Writing Input File for RSMinerve
The input file format for RSMinerve is described in Garcia Hernandez et al. (2020), page 136.
Q <- Qimb_m3s_sub |>
mutate(Qimb_m3s = round(Qimb_m3s, digits = 7))
temp_wide <- Q |>
pivot_wider(names_from = name_2, values_from = Qimb_m3s) |>
rename(Station = date)
datechar <- posixct2rsminerveChar(temp_wide$Station)$value
datechar <- gsub(" 01:00:00", " 00:00:00", datechar)
datechar <- gsub(" 02:00:00", " 00:00:00", datechar)
output <- rbind(colnames(temp_wide),
c("X", "1", "1", "1", "1"), # Random coordinates, not relevant
c("Y", "2", "2", "2", "3"),
c("Z", "3", "3", "3", "3"),
c("Sensor", "Q", "Q", "Q", "Q"),
c("Category", "Flow", "Flow", "Flow", "Flow"),
c("Unit", "m3/s", "m3/s", "m3/s", "m3/s"),
c("Interpolation", "Linear", "Linear", "Linear",
"Linear"),
cbind(datechar,
as.character(temp_wide$Dzhaldzhur_Subbasin),
as.character(temp_wide$Ulak_Subbasin),
as.character(temp_wide$Atbaschy_Midstream_Subbasin),
as.character(temp_wide$Atbaschy_Downstream_Subbasin)))Finally, the prepared data is written to a csv file which can be read into RSMinerve.
writefilename <- paste0(data_path, "RSM_demo_glacier_source.csv")
write.table(output, file = writefilename, col.names = FALSE,
row.names = FALSE, append = FALSE, quote = FALSE,
sep = ",", dec = ".")12.2 Implementation in RS Minerve
In RSMinerve, a source object is placed for each sub-basin containing glaciers. The RSMinerve input file written above is then read into the database and connected to the appropriate source objects in the model window. The RSMinerve model including glacier discharge can now be run.
12.3 References
Erasov, N. V. 1968. “Method for Determining of Volume of Mountain Glaciers.” MGI, no. 14: 307–8.
Farinotti, Daniel, Matthias Huss, Johannes J. Fürst, Johannes Landmann, Horst Machguth, Fabien Maussion, and Ankur Pandit. 2019. “A Consensus Estimate for the Ice Thickness Distribution of All Glaciers on Earth.” Nature Geoscience 12 (3): 168–73. https://doi.org/10.1038/s41561-019-0300-3.
Garcia Hernandez, J., A. Foehn, J. Fluixa-Sanmartin, B. Roquier, T. Brauchli, J. Paredes Arquiola, and De Cesare G. 2020. “RS MINERVE - Technical Manual, V2.25.” ISSN 2673-2661. Switzerland: Ed. CREALP.
Hock, Regine. 2003. “Temperature index melt modelling in mountain areas.” Journal of Hydrology 282 (1-4): 104–15. https://doi.org/10.1016/s0022-1694(03)00257-9.
Hugonnet, Romain, Robert McNabb, Etienne Berthier, Brian Menounos, Christopher Nuth, Luc Girod, Daniel Farinotti, et al. 2021. “Accelerated Global Glacier Mass Loss in the Early Twenty-First Century.” Nature 592 (7856): 726–31. https://doi.org/10.1038/s41586-021-03436-z.
Karger, Dirk Nikolaus, Olaf Conrad, Jürgen Böhner, Tobias Kawohl, Holger Kreft, Rodrigo Wilber Soria-Auza, Niklaus E. Zimmermann, H. Peter Linder, and Michael Kessler. 2017. “Climatologies at high resolution for the earth’s land surface areas.” Scientific Data 4 (1): 170122. https://doi.org/10.1038/sdata.2017.122.
Miles, Evan, Michael McCarthy, Amaury Dehecq, Marin Kneib, Stefan Fugger, and Francesca Pellicciotti. 2021. “Health and Sustainability of Glaciers in High Mountain Asia.” Nature Communications 12 (2868): 10. https://doi.org/https://doi.org/10.1038/s41467-021-23073-4.